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Abstract 

An iterative algorithm is presented for solving the RPA equations of linear 
response. The method optimally computes the energy-weighted moments of 
the strength function, allowing one to match the computational effort to the 
intrinsic accuracy of the basic mean- field approximation, avoiding the problem 
of solving very large matrices. For local interactions, the computational effort 
for the method scales with the number of particles N p as O(Np). 

I. INTRODUCTION 

In a number of branches of physics, mean field theory gives a remarkably effective approx- 
imation to the ground state. Similarly, for the response of the system to small perturbations, 
time-dependent mean field theory is a useful extension. This is the experience in nuclear 
physics [|I],f2[ , atomic and molecular physics [||-|8j and condensed matter physics PJITl . There 



are of course intrinsic limitations to these approximations, but equally pressing is the large 
computational resources required for calculations of systems of interest. This is our motiva- 
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tion to look for algorithms that better match the computational effort to the intrinsic limits 
of the approximation. 



We take our inspiration from the Lanczos algorithm [11], which is best known in many 



body physics for extracting low-lying eigenstates of very large Hamiltonian matrices ||12| , |T3 
When dealing with large spaces, the computational question often comes down to the number 
of times the Hamiltonian operates on a state vector. Depending on the Hamiltonian and the 
starting vector, the Lanczos algorithm is able to extract an accurate ground state vector in 
a basis of 10 5-6 states with a few hundred Hamiltonian operations. The algorithm may be 
viewed as a numerically stable technique |14| to compute moments of the Hamiltonian with 
respect to some initial state that is, ^ = (^o| H k |\l>o) • F° r large k, /i^ is dominated by 



the extremal eigenvalues [15], which are thus available for recovery. 



The Lanczos algorithm has also been applied to many other topics in atomic, molecular, 



solid state, and nuclear physics, including computation of the S-matrix [16], time-evolution of 



wave packets |T7|], level densities |18[, and the continued- fraction expansion of the resolvant 



or Green's function Particularly relevant to us is the application to strength functions. 
The strength function S for an operator Q on a state i is defined 

s(E) = Y< 6 ( E - E f + E >)\(f\Q\*)\ 2 - (!) 

/ 

A powerful technique to calculate the strength function, successfully applied to the nu- 



clear shell model ||T3|pO| , uses the Lanczos algorithm with a starting vector \^ ) = 
Q \i) / {i\Q 2 \i) l l 2 ■ The Lanczos algorithm implicitly computes the moments 

M k = J dE(E - E t ) k S(E) 

of the strength function. After a few tens of iterations one can accurately reconstruct the 
distribution of the exact strength function. 

In many cases, however, the matrix elements of the operator Q are sensitive to correla- 
tions in the ground state, and then the size of the wave function basis in the straightforward 
Hamiltonian approach becomes problematic. In this situation, the time-dependent mean- 
field theory offers a reasonable compromise. The small amplitude theory, the RPA or linear 



response, can be cast into a matrix form in a particle-hole basis. However, the RPA matrix 
is not symmetric as required by the Lanczos algorithm. The matrix equation is commonly 
written as 



UJ 



(2) 



A B \ x 

-B -AJ \yj \y, 
where A and B are particle-hole Hamiltonian matrices, u is the eigenfrequency, and x and 

y are the vectors of positive- and negative-frequency particle-hole amplitudes, respectively. 

An important property of the RPA equation is that eigenvectors come in conjugate pairs: in 

equation ([2J) (y, x) is also an eigenvector with eigenfrequency — to. For the linear response, 

the matrix element between the RPA ground state |0) and an excited state \u) may be 

expressed as 



{u\Q\0) = (q,q) 



x 



V, 



(3) 



where q is the vector of particle-hole matrix elements and the vector (x, y) is normalized as 
1 = x ■ x — y ■ y. 

There are a number of ways to introduce a Lanczos-type algorithm for the RPA matrix. 

The method we describe here has the advantages that it preserves the form eq.(§) of the 

RPA matrix and it produces strength functions that respect sum rules. We seek a new basis 

of vectors \Zj) : = (Aj,Fj) where the matrices of column vectors U := (Xi, X 2 , X 3 , . . .) and 

V := (Yi, Y 2 , ^3) ■ ■ ■) transform the RPA matrix as 

U T -V T \ / A B WU V\ /A' 

-V T U T ) V-B -A / \ V U;~ V-B' 
where the transformed matrices A' and B' are now tridiagonal: 

/ e% a\ \ / d\ bi 



B' 

-A' 



(4) 



A' 



ai e 2 a 2 
a 2 e 3 



V 



B' 



bi d 2 b 2 
b 2 d 3 



V 



(5) 



The Lanczos basis vectors and matrix elements are generated iteratively as follows. Sup- 
pose we have the vectors \Zx), \Z n ) already computed, together with the transformed 



matrix up to e n _i, d n -i, a n -\ and 6 n -i- The iteration starts by applying the RPA matrix in 
eqn. (0) to the vector \Z n ), 

(X t \ ( AX n + BY n \ 
\Zt) = k = . J (6) 

Vy t y \-Bx n -AY n ) 

The diagonal elements e n and d n are now easily computed: 

&n — Xf X n — YfY n 

d n = X t -Y n — Y t ■ X n . (7) 

We next project out \Z£), the component of \Z t ) that is orthogonal to the space \Z n ). 
This can be done conveniently by using the matrix elements in (|^) that have already been 
calculated, 

(X t \ ( Xt — e n X n + d n Y n — a n -iX n _i + b 
= . . (8) 

/ V d n X n -\- e n Y n 

The norm of the vector is then computed as 

A/- = x;-x;-f/-f/ (9) 

The norm can be negative, and the definition of the new vector |Z n +i) depends on the sign. 
In fact, because we are actually doing 6/oc/i;-Lanczos, implicitly operating not only on the 
vector(X, Y) but also its RPA conjugate (Y, X) simultaneously, there is a degree of freedom, 
corresponding to a hyperbolic rotation, in choosing the new vector. The simplest choice for 
the vectors and corresponding RPA matrix elements is 
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'XI 



VAT \ Y! 



r n+1 ) = -y=\ _ , a n+1 = VAf, b n+1 = 0; Af > (10) 



and 



\Z n+1 ) = -j=j II, a n+l = 0, b n+1 = V^S, J\f < (11) 

This completes the iteration cycle. 

In analogy with the application to strength functions in the nuclear shell model, we start 
with the vector given by 



(12) 

With such a starting vector the algorithm manifestly preserves the energy-weighted sum 
rules: 

M k = J2^(^\Q\0}\ k odd. (13) 

V 

Using the eigenvector representation of the RPA matrix, one can show 

M " = \&® ( \ ^) (^) ?A;0dd - (M) 
With our method the n-th iterate respects the odd-A; sum rules for k < 2n — 1. 

We now illustrate the method with a very simple model, a collective particle-hole inter- 
action fragmented by single-particle energies. We consider states i = 1, ...,JV with matrix 
elements A^- = + nqiqj and = Kq^qj. Here e represents the energy spacing of the 
particle-hole configurations, k is the strength of the collective coupling to the field Q, and 
the components of the vector q { oc i(N — i) x r, where the r are Gaussian distributed ran- 
dom amplitudes, and normalize \q\ 2 = 1. The factor i(N — i) weights the collective response 
towards the middle of the excitation spectrum. The parameter k should be positive for a 
repulsive collective interaction such as the Coulomb that generates plasmons. 

In Fig. 1 we show the strength function for such an RPA matrix in a space of 500 
states, with parameter values given in the caption. The parameters were chosen to obtain 
moderate collectivity, with a strong but broadly-fragmented collective excitation distributed 
over the spectrum. Fig. 1 also displays the n=3, 10, and 50 approximants to the strength 
function, where n is the number of Lanczos vectors \Zj), or, equivalently, the number of 
multiplications with the RPA matrix. One sees that with a handful of states, one state 
closely approximates the collective excitation and the others distribute themselves over the 
remaining spectrum. A better way to see the convergence of the strength function is to 
plot its integral, = I]^6(w — cj^)(^|(5|0) 2 . This is shown in Fig. 2 for n=3 and 10. 
After 50 iterations the integral of the strength function is virtually indistinguishable from 
the exact solution. 
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We mention that the algorithm does not explicitly preserve the total strength Mq. If 
there were no correlations in the ground state, that is, if the vectors Yj all vanished, then 
the total strength would be \q\ 2 = 1. The non-trivial deviations from 1 in our examples 
are related to the amount of correlations in the ground state. This is illustrated in Fig. 3, 
which is the integrated strength function for a model identical to that in Figs. 1,2 except 
that the collective interaction is attractive rather than repulsive. Here the total strength 
is about 3.7, i.e. quite different from 1. Fig. 3 illustrates how with n =3 and 10 the total 
approximate strength converges rapidly to the exact value. (In the repulsive model of Fig. 2 
the total strength had already converged by n = 3.) Although we cannot prove this rapid 
convergence in all cases, it seems likely in light of the strong constraints imposed by the 
odd-k sum rules. 

We anticipate that the algorithm will be particularly useful in problems which require a 
single-particle dimensionality of the order of tens or hundreds of thousands, but which allow 
a sparse matrix approximation for the Hamiltonian, such as the local density approxima- 
tion. This applies to molecular and condensed matter physics modeled with the Kohn-Sham 
equations, and to nuclear physics for excitations in deformed nuclei 0. With the LDA 
Hamiltonian, an efficient particle-hole representation can be constructed from the orbital 
representation of holes and the coordinate-space representation of particles fl21 |. The com- 



putational difficulty for the basic matrix-vector multiplication then scales as the number of 
particles N p and the dimensionality of the single-particle space M as MN£ ~ iV?. Only 
a fixed number of these operations, of the order of ten, are needed to obtain the strength 
function to the accuracy of the fundamental mean field approximation. Thus the overall 
scaling of the method is O(A^). 

This study arose in the program at the Institute for Nuclear Theory, "Numerical methods 
for strongly interacting quantum systems", and we wish to thank J. Carlson and R. Wiringa 
for providing that forum. G.B. also thanks K. Yabana for many discussions. Financial 
support was provided by the INT under Department of Energy Grant FG06-90ER40561 
and by Department of Energy Grant DE-FG02-96ER40985. 



FIGURES 

FIG. 1. Strength function for the model described in the text with k = 10 and e = 0.1 (in 
arbitrary units) for 500 states, and the Lanczos approximants for 5, 10, and 50 Lanczos vectors. 
The scales for the abscissae are different because the strength is fragmented over a different number 
of states. 

FIG. 2. Integrated strength function for the model described in Fig. 1. For 50 Lanczos vectors 
the integrated strength is virtually indistinguishable on this graph from the full calculation. 

FIG. 3. The same as figure 2, except with the collective interaction is attractive, k = — 10. 
Notice that the total strength is not constrained, as described in the text, but has converged by 
the 10th iteration. 
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